/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.linear;

import bd.org.apache.commons.math.exception.*;
import bd.org.apache.commons.math.exception.util.LocalizedFormats;
import bd.org.apache.commons.math.util.FastMath;
import bd.org.apache.commons.math.util.Precision;

import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;

/**
 * A collection of static methods that operate on or return matrices.
 *
 * @version $Id: MatrixUtils.java 1244107 2012-02-14 16:17:55Z erans $
 */
public class MatrixUtils {

    /**
     * Private constructor.
     */
    private MatrixUtils() {

        super();
    }

    /**
     * Returns a {@link RealMatrix} with specified dimensions.
     * <p>The type of matrix returned depends on the dimension. Below
     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
     * square matrix) which can be stored in a 32kB array, a {@link
     * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
     * BlockRealMatrix} instance is built.</p>
     * <p>The matrix elements are all set to 0.0.</p>
     *
     * @param rows number of rows of the matrix
     * @param columns number of columns of the matrix
     * @return RealMatrix with specified dimensions
     * @see #createRealMatrix(double[][])
     */
    public static RealMatrix createRealMatrix(final int rows, final int columns) {

        return (rows * columns <= 4096)
               ? new Array2DRowRealMatrix(rows, columns)
               : new BlockRealMatrix(rows, columns);
    }


    /**
     * Returns a {@link RealMatrix} whose entries are the the values in the
     * the input array.
     * <p>The type of matrix returned depends on the dimension. Below
     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
     * square matrix) which can be stored in a 32kB array, a {@link
     * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
     * BlockRealMatrix} instance is built.</p>
     * <p>The input array is copied, not referenced.</p>
     *
     * @param data input array
     * @return RealMatrix containing the values of the array
     * is {@code null}.
     * @see #createRealMatrix(int, int)
     */
    public static RealMatrix createRealMatrix(double[][] data) {

        if ((data == null) || (data[0] == null)) {
            throw new NullArgumentException();
        }

        return (data.length * data[0].length <= 4096)
               ? new Array2DRowRealMatrix(data)
               : new BlockRealMatrix(data);
    }


    /**
     * Returns <code>dimension x dimension</code> identity matrix.
     *
     * @param dimension dimension of identity matrix to generate
     * @return identity matrix
     * @since 1.1
     */
    public static RealMatrix createRealIdentityMatrix(int dimension) {

        final RealMatrix m = createRealMatrix(dimension, dimension);

        for (int i = 0; i < dimension; ++i) {

            m.setEntry(i, i, 1.0);
        }

        return m;
    }


    /**
     * Returns a diagonal matrix with specified elements.
     *
     * @param diagonal diagonal elements of the matrix (the array elements
     * will be copied)
     * @return diagonal matrix
     * @since 2.0
     */
    public static RealMatrix createRealDiagonalMatrix(final double[] diagonal) {

        final RealMatrix m = createRealMatrix(diagonal.length, diagonal.length);

        for (int i = 0; i < diagonal.length; ++i) {

            m.setEntry(i, i, diagonal[i]);
        }

        return m;
    }


    /**
     * Creates a {@link RealVector} using the data from the input array.
     *
     * @param data the input data
     * @return a data.length RealVector
     */
    public static RealVector createRealVector(double[] data) {

        if (data == null) {
            throw new NullArgumentException();
        }

        return new ArrayRealVector(data, true);
    }


    /**
     * Create a row {@link RealMatrix} using the data from the input
     * array.
     *
     * @param rowData the input row data
     * @return a 1 x rowData.length RealMatrix
     */
    public static RealMatrix createRowRealMatrix(double[] rowData) {

        if (rowData == null) {
            throw new NullArgumentException();
        }

        final int        nCols = rowData.length;
        final RealMatrix m     = createRealMatrix(1, nCols);

        for (int i = 0; i < nCols; ++i) {

            m.setEntry(0, i, rowData[i]);
        }

        return m;
    }


    /**
     * Creates a column {@link RealMatrix} using the data from the input
     * array.
     *
     * @param columnData the input column data
     * @return a columnData x 1 RealMatrix
     */
    public static RealMatrix createColumnRealMatrix(double[] columnData) {

        if (columnData == null) {
            throw new NullArgumentException();
        }

        final int        nRows = columnData.length;
        final RealMatrix m     = createRealMatrix(nRows, 1);

        for (int i = 0; i < nRows; ++i) {

            m.setEntry(i, 0, columnData[i]);
        }

        return m;
    }


    /**
     * Check if matrix indices are valid.
     *
     * @param m Matrix.
     * @param row Row index to check.
     * @param column Column index to check.
     */
    public static void checkMatrixIndex(final AnyMatrix m, final int row, final int column) {

        checkRowIndex(m, row);
        checkColumnIndex(m, column);
    }


    /**
     * Check if a row index is valid.
     *
     * @param m Matrix.
     * @param row Row index to check.
     */
    public static void checkRowIndex(final AnyMatrix m, final int row) {

        if ((row < 0) || (row >= m.getRowDimension())) {

            throw new OutOfRangeException(LocalizedFormats.ROW_INDEX, row, 0,
                                          m.getRowDimension() - 1);
        }
    }


    /**
     * Check if a column index is valid.
     *
     * @param m Matrix.
     * @param column Column index to check.
     */
    public static void checkColumnIndex(final AnyMatrix m, final int column) {

        if ((column < 0) || (column >= m.getColumnDimension())) {

            throw new OutOfRangeException(LocalizedFormats.COLUMN_INDEX, column, 0,
                                          m.getColumnDimension() - 1);
        }
    }


    /**
     * Check if submatrix ranges indices are valid.
     * Rows and columns are indicated counting from 0 to {@code n - 1}.
     *
     * @param m Matrix.
     * @param startRow Initial row index.
     * @param endRow Final row index.
     * @param startColumn Initial column index.
     * @param endColumn Final column index.
     * {@code endColumn < startColumn}.
     */
    public static void checkSubMatrixIndex(final AnyMatrix m, final int startRow, final int endRow,
            final int startColumn, final int endColumn) {

        checkRowIndex(m, startRow);
        checkRowIndex(m, endRow);

        if (endRow < startRow) {

            throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_ROW_AFTER_FINAL_ROW,
                    endRow, startRow, false);
        }

        checkColumnIndex(m, startColumn);
        checkColumnIndex(m, endColumn);

        if (endColumn < startColumn) {

            throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
                    endColumn, startColumn, false);
        }

    }


    /**
     * Check if submatrix ranges indices are valid.
     * Rows and columns are indicated counting from 0 to n-1.
     *
     * @param m Matrix.
     * @param selectedRows Array of row indices.
     * @param selectedColumns Array of column indices.
     * {@code selectedColumns} are {@code null}.
     */
    public static void checkSubMatrixIndex(final AnyMatrix m, final int[] selectedRows,
            final int[] selectedColumns) {

        if (selectedRows == null) {
            throw new NullArgumentException();
        }

        if (selectedColumns == null) {
            throw new NullArgumentException();
        }

        if (selectedRows.length == 0) {
            throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
        }

        if (selectedColumns.length == 0) {
            throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_COLUMN_INDEX_ARRAY);
        }

        for (final int row : selectedRows) {

            checkRowIndex(m, row);
        }

        for (final int column : selectedColumns) {

            checkColumnIndex(m, column);
        }
    }


    /**
     * Check if matrices are addition compatible.
     *
     * @param left Left hand side matrix.
     * @param right Right hand side matrix.
     */
    public static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right) {

        if ((left.getRowDimension() != right.getRowDimension())
                || (left.getColumnDimension() != right.getColumnDimension())) {

            throw new MatrixDimensionMismatchException(left.getRowDimension(),
                    left.getColumnDimension(), right.getRowDimension(), right.getColumnDimension());
        }
    }


    /**
     * Check if matrices are subtraction compatible
     *
     * @param left Left hand side matrix.
     * @param right Right hand side matrix.
     */
    public static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right) {

        if ((left.getRowDimension() != right.getRowDimension())
                || (left.getColumnDimension() != right.getColumnDimension())) {

            throw new MatrixDimensionMismatchException(left.getRowDimension(),
                    left.getColumnDimension(), right.getRowDimension(), right.getColumnDimension());
        }
    }


    /**
     * Check if matrices are multiplication compatible
     *
     * @param left Left hand side matrix.
     * @param right Right hand side matrix.
     */
    public static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right) {

        if (left.getColumnDimension() != right.getRowDimension()) {

            throw new DimensionMismatchException(left.getColumnDimension(),
                    right.getRowDimension());
        }
    }


    /**
     * Serialize a {@link RealVector}.
     * <p>
     * This method is intended to be called from within a private
     * <code>writeObject</code> method (after a call to
     * <code>oos.defaultWriteObject()</code>) in a class that has a
     * {@link RealVector} field, which should be declared <code>transient</code>.
     * This way, the default handling does not serialize the vector (the {@link
     * RealVector} interface is not serializable by default) but this method does
     * serialize it specifically.
     * </p>
     * <p>
     * The following example shows how a simple class with a name and a real vector
     * should be written:
     * <pre><code>
     * public class NamedVector implements Serializable {
     *     private final String name;
     *     private final transient RealVector coefficients;
     *     // omitted constructors, getters ...
     *     private void writeObject(ObjectOutputStream oos) throws IOException {
     *         oos.defaultWriteObject();  // takes care of name field
     *         MatrixUtils.serializeRealVector(coefficients, oos);
     *     }
     *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
     *         ois.defaultReadObject();  // takes care of name field
     *         MatrixUtils.deserializeRealVector(this, "coefficients", ois);
     *     }
     * }
     * </code></pre>
     * </p>
     *
     * @param vector real vector to serialize
     * @param oos stream where the real vector should be written
     * @throws IOException if object cannot be written to stream
     * @see #deserializeRealVector(Object, String, ObjectInputStream)
     */
    public static void serializeRealVector(final RealVector vector, final ObjectOutputStream oos)
            throws IOException {

        final int n = vector.getDimension();

        oos.writeInt(n);

        for (int i = 0; i < n; ++i) {

            oos.writeDouble(vector.getEntry(i));
        }
    }


    /**
     * Deserialize  a {@link RealVector} field in a class.
     * <p>
     * This method is intended to be called from within a private
     * <code>readObject</code> method (after a call to
     * <code>ois.defaultReadObject()</code>) in a class that has a
     * {@link RealVector} field, which should be declared <code>transient</code>.
     * This way, the default handling does not deserialize the vector (the {@link
     * RealVector} interface is not serializable by default) but this method does
     * deserialize it specifically.
     * </p>
     *
     * @param instance instance in which the field must be set up
     * @param fieldName name of the field within the class (may be private and final)
     * @param ois stream from which the real vector should be read
     * @throws ClassNotFoundException if a class in the stream cannot be found
     * @throws IOException if object cannot be read from the stream
     * @see #serializeRealVector(RealVector, ObjectOutputStream)
     */
    public static void deserializeRealVector(final Object instance, final String fieldName,
            final ObjectInputStream ois)
            throws ClassNotFoundException, IOException {

        try {

            // read the vector data
            final int      n    = ois.readInt();
            final double[] data = new double[n];

            for (int i = 0; i < n; ++i) {

                data[i] = ois.readDouble();
            }

            // create the instance
            final RealVector vector = new ArrayRealVector(data, false);

            // set up the field
            final java.lang.reflect.Field f = instance.getClass().getDeclaredField(fieldName);

            f.setAccessible(true);
            f.set(instance, vector);

        } catch (NoSuchFieldException nsfe) {

            IOException ioe = new IOException();

            ioe.initCause(nsfe);

            throw ioe;

        } catch (IllegalAccessException iae) {

            IOException ioe = new IOException();

            ioe.initCause(iae);

            throw ioe;
        }

    }


    /**
     * Serialize a {@link RealMatrix}.
     * <p>
     * This method is intended to be called from within a private
     * <code>writeObject</code> method (after a call to
     * <code>oos.defaultWriteObject()</code>) in a class that has a
     * {@link RealMatrix} field, which should be declared <code>transient</code>.
     * This way, the default handling does not serialize the matrix (the {@link
     * RealMatrix} interface is not serializable by default) but this method does
     * serialize it specifically.
     * </p>
     * <p>
     * The following example shows how a simple class with a name and a real matrix
     * should be written:
     * <pre><code>
     * public class NamedMatrix implements Serializable {
     *     private final String name;
     *     private final transient RealMatrix coefficients;
     *     // omitted constructors, getters ...
     *     private void writeObject(ObjectOutputStream oos) throws IOException {
     *         oos.defaultWriteObject();  // takes care of name field
     *         MatrixUtils.serializeRealMatrix(coefficients, oos);
     *     }
     *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
     *         ois.defaultReadObject();  // takes care of name field
     *         MatrixUtils.deserializeRealMatrix(this, "coefficients", ois);
     *     }
     * }
     * </code></pre>
     * </p>
     *
     * @param matrix real matrix to serialize
     * @param oos stream where the real matrix should be written
     * @throws IOException if object cannot be written to stream
     * @see #deserializeRealMatrix(Object, String, ObjectInputStream)
     */
    public static void serializeRealMatrix(final RealMatrix matrix, final ObjectOutputStream oos)
            throws IOException {

        final int n = matrix.getRowDimension();
        final int m = matrix.getColumnDimension();

        oos.writeInt(n);
        oos.writeInt(m);

        for (int i = 0; i < n; ++i) {

            for (int j = 0; j < m; ++j) {

                oos.writeDouble(matrix.getEntry(i, j));
            }
        }
    }


    /**
     * Deserialize  a {@link RealMatrix} field in a class.
     * <p>
     * This method is intended to be called from within a private
     * <code>readObject</code> method (after a call to
     * <code>ois.defaultReadObject()</code>) in a class that has a
     * {@link RealMatrix} field, which should be declared <code>transient</code>.
     * This way, the default handling does not deserialize the matrix (the {@link
     * RealMatrix} interface is not serializable by default) but this method does
     * deserialize it specifically.
     * </p>
     *
     * @param instance instance in which the field must be set up
     * @param fieldName name of the field within the class (may be private and final)
     * @param ois stream from which the real matrix should be read
     * @throws ClassNotFoundException if a class in the stream cannot be found
     * @throws IOException if object cannot be read from the stream
     * @see #serializeRealMatrix(RealMatrix, ObjectOutputStream)
     */
    public static void deserializeRealMatrix(final Object instance, final String fieldName,
            final ObjectInputStream ois)
            throws ClassNotFoundException, IOException {

        try {

            // read the matrix data
            final int        n    = ois.readInt();
            final int        m    = ois.readInt();
            final double[][] data = new double[n][m];

            for (int i = 0; i < n; ++i) {

                final double[] dataI = data[i];

                for (int j = 0; j < m; ++j) {

                    dataI[j] = ois.readDouble();
                }
            }

            // create the instance
            final RealMatrix matrix = new Array2DRowRealMatrix(data, false);

            // set up the field
            final java.lang.reflect.Field f = instance.getClass().getDeclaredField(fieldName);

            f.setAccessible(true);
            f.set(instance, matrix);

        } catch (NoSuchFieldException nsfe) {

            IOException ioe = new IOException();

            ioe.initCause(nsfe);

            throw ioe;

        } catch (IllegalAccessException iae) {

            IOException ioe = new IOException();

            ioe.initCause(iae);

            throw ioe;
        }
    }


    /**
     * Solve  a  system of composed of a Lower Triangular Matrix
     * {@link RealMatrix}.
     * <p>
     * This method is called to solve systems of equations which are
     * of the lower triangular form. The matrix {@link RealMatrix}
     * is assumed, though not checked, to be in lower triangular form.
     * The vector {@link RealVector} is overwritten with the solution.
     * The matrix is checked that it is square and its dimensions match
     * the length of the vector.
     * </p>
     *
     * @param rm RealMatrix which is lower triangular
     * @param b RealVector this is overwritten
     */
    public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b) {

        if ((rm == null) || (b == null) || (rm.getRowDimension() != b.getDimension())) {

            throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
                    (rm == null)
                    ? 0
                    : rm.getRowDimension(), (b == null)
                                            ? 0
                                            : b.getDimension());
        }

        if (rm.getColumnDimension() != rm.getRowDimension()) {

            throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
                    rm.getRowDimension(), rm.getRowDimension(), rm.getRowDimension(),
                    rm.getColumnDimension());
        }

        int rows = rm.getRowDimension();

        for (int i = 0; i < rows; i++) {

            double diag = rm.getEntry(i, i);

            if (FastMath.abs(diag) < Precision.SAFE_MIN) {
                throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
            }

            double bi = b.getEntry(i) / diag;

            b.setEntry(i, bi);

            for (int j = i + 1; j < rows; j++) {

                b.setEntry(j, b.getEntry(j) - bi * rm.getEntry(j, i));
            }
        }
    }


    /**
     * Solver a  system composed  of an Upper Triangular Matrix
     * {@link RealMatrix}.
     * <p>
     * This method is called to solve systems of equations which are
     * of the lower triangular form. The matrix {@link RealMatrix}
     * is assumed, though not checked, to be in upper triangular form.
     * The vector {@link RealVector} is overwritten with the solution.
     * The matrix is checked that it is square and its dimensions match
     * the length of the vector.
     * </p>
     *
     * @param rm RealMatrix which is upper triangular
     * @param b RealVector this is overwritten
     */
    public static void solveUpperTriangularSystem(RealMatrix rm, RealVector b) {

        if ((rm == null) || (b == null) || (rm.getRowDimension() != b.getDimension())) {

            throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
                    (rm == null)
                    ? 0
                    : rm.getRowDimension(), (b == null)
                                            ? 0
                                            : b.getDimension());
        }

        if (rm.getColumnDimension() != rm.getRowDimension()) {

            throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
                    rm.getRowDimension(), rm.getRowDimension(), rm.getRowDimension(),
                    rm.getColumnDimension());
        }

        int rows = rm.getRowDimension();

        for (int i = rows - 1; i > -1; i--) {

            double diag = rm.getEntry(i, i);

            if (FastMath.abs(diag) < Precision.SAFE_MIN) {
                throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
            }

            double bi = b.getEntry(i) / diag;

            b.setEntry(i, bi);

            for (int j = i - 1; j > -1; j--) {

                b.setEntry(j, b.getEntry(j) - bi * rm.getEntry(j, i));
            }
        }
    }
}
